## Admissions-as-corrections reduce support for partisan misperceptions 
#and intended partisan media consumption
##(Ms. No. 240015)--The Journal of Politics##
##Replication code for JOP Dataverse##
##Prepared by J.J. Fahey, 3/20/2025##

##The following libraries are necessary for analysis and data viz##
library(rio)
library(foreign)
library(tidyverse)
library(tidyr)
library(dplyr)
library(estimatr)
library(lavaan)
library(modelsummary)
library(ggplot2)
library(broom)
library(TOSTER)

#Set working directory and import data


aac.data <- import("JOP_AAC.csv")

####Experiment order and treatment conditions####
##Note: Sandy Hook is considered the "first," experiment
#Fox News is "second" 

##Experiment on Sandy Hook## 
##Converting qualtrics randomization IDs to numeric##
aac.data <- aac.data %>%
  mutate(treatment.jones = case_when (
    aac.data$FL_61_DO == "FL_62" ~ 0, #control
    aac.data$FL_61_DO  == "FL_64" ~ 1, #admissions
    aac.data$FL_61_DO == "FL_65" ~ 2)) #fact check

##Experiment on Fox News## 
aac.data <- aac.data %>%
  mutate(treatment.fox = case_when (
    aac.data$FL_22_DO == "FL_23" ~ 0, #control
    aac.data$FL_22_DO == "FL_24" ~ 1, #admissions
    aac.data$FL_22_DO== "FL_63" ~ 2)) #fact check

####Data quality checks####

aac.data <- aac.data %>%
  mutate (complete = case_when (
    aac.data$Progress == 100 ~ 1,
    aac.data$Progress < 100 ~ 0
  ))

##Note: to calculate differential attrition later (Line 879), do not run 
##lines 55-80 as they drop out of sample/inattentive respondents

aac.data <- aac.data[aac.data$complete == 1, ]

##Dropping respondents who indicate they are democrats##
aac.data <- aac.data[aac.data$ind1 == 0 | aac.data$ind1 == 2 | aac.data$pid == 1,]

##Removing inattentive respondents##
##Note: Attn Check 2 removed from pilot for being insufficiently discriminatory##

##Attn Check 1##
aac.data <- aac.data %>%
  mutate(attncheck1= case_when (
    aac.data$attn1 == 4 ~ 1, 
    TRUE ~ 0))

##Attn Check 3##
aac.data <- aac.data %>%
  mutate(attncheck3= case_when (
    aac.data$attn3 == 1 ~ 1, 
    TRUE ~ 0 ))

##Dropping respondents who failed more than one attncheck 
aac.data <- aac.data %>%
  mutate(attntotal = attncheck1 + attncheck3)

aac.data <- aac.data[aac.data$attntotal >= 1, ]

####Pre-Analysis Demographic Manipulation####
##Change age to numeric variable

aac.data$age <- as.numeric(aac.data$age)
is.numeric(aac.data$age)

table(aac.data$age)

##Four point PID scale (because Dems are dropped)##
aac.data <- aac.data %>%
  mutate(pid.scale = case_when(
    aac.data$rep1 == 1~1, 
    aac.data$rep1 == 0~.666, 
    aac.data$ind1 == 0 ~ .333,
    aac.data$ind1 == 2 ~ 0))

##White dummy##
##Note: only considered white if they select *only* white 
aac.data <- aac.data %>%
  mutate(white = case_when(
    aac.data$race == 0~1,
    TRUE ~ 0))

table(aac.data$white)

##Hispanic/Latino dummy##
aac.data <- aac.data %>%
  mutate(hispanic = case_when(
    aac.data$hispanic == 1~1,
    TRUE ~ 0))

table(aac.data$hispanic)

##College dummmy## 

aac.data <- aac.data %>%
  mutate(college = case_when (
    aac.data$educ <= 4 ~ 0, 
    aac.data$educ >= 5 ~ 1
  ))

##Male Dummy## 
aac.data <- aac.data %>%
  mutate(male = case_when (
    aac.data$gender == 0 ~1, 
    aac.data$gender == 1 ~0, 
    aac.data$gender == 2 ~0))

##Media trust## 
aac.data <- aac.data %>%
  mutate(media.trust = case_when (
    aac.data$mediatrust== 0 ~ 0,
    aac.data$mediatrust == 1 ~ .25,
    aac.data$mediatrust == 2 ~ .5,
    aac.data$mediatrust == 3 ~ .75,
    aac.data$mediatrust== 4 ~ 1))

##Fox Viewer##
##Included as Fox viewer if they select Fox at all
aac.data <- aac.data %>%
  mutate(fox = case_when(
    str_detect(tvstations, "1") ~ 1,
    TRUE ~ 0
  ))

table(aac.data$fox)

##Jones Viewer## 
##Included as Jones viewer if they select Jones at all
aac.data <- aac.data %>%
  mutate(jones = case_when(
    str_detect(alternatepods, "1") ~ 1,
    TRUE ~ 0))

table(aac.data$jones)

####Creating dependent variables#### 
##Note: H1a & H1b variables measured both pre and post...
#...and included as controls in models 

##Conspiracy variables (H1a and H1b)## 
##Trump won (pre)##
aac.data <- aac.data %>%
  mutate(electionwinner.pre = case_when (
    aac.data$electionwinnerpre == 1 ~ 1, 
    aac.data$electionwinnerpre == 2 ~ 0, 
    aac.data$electionwinnerpre == 0 ~ 0))

table(aac.data$electionwinner.pre)

##Dominion changed votes (pre)## 
aac.data <- aac.data %>%
  mutate(voteschanged.pre = case_when (
    aac.data$voteschangedpre == 1 ~ 1, 
    aac.data$voteschangedpre == 2 ~ 0, 
    aac.data$voteschangedpre == 0 ~ 0))

table(aac.data$voteschanged.pre)

##Immigrants Voted (pre)## 
aac.data <- aac.data %>%
  mutate(illegalvotes.pre = case_when (
    aac.data$illegalvotespre == 1 ~ 1, 
    aac.data$illegalvotespre == 2 ~ 0, 
    aac.data$illegalvotespre == 0 ~ 0))

table(aac.data$illegalvotes.pre)

##venezuela (pre)## 
aac.data <- aac.data %>%
  mutate(venezuela.pre= case_when (
    aac.data$venezuelapre == 0 ~ 1, 
    aac.data$venezuelapre == 2 ~ 0, 
    aac.data$venezuelapre == 1 ~ 0))

table(aac.data$venezuela.pre)

##uvalde (pre)## 
aac.data <- aac.data %>%
  mutate(uvalde.pre= case_when (
    aac.data$uvaldepre == 0 ~ 0, 
    aac.data$uvaldepre == 2 ~ 0, 
    aac.data$uvaldepre == 1 ~ 1))

##sandy hook (pre)## 
aac.data <- aac.data %>%
  mutate(sandyhook.pre= case_when (
    aac.data$sandyhookpre == 1 ~ 1, 
    aac.data$sandyhookpre == 2 ~ 0, 
    aac.data$sandyhookpre == 0 ~ 0))

##Post stimuli variables##

##Trump won (post)##
aac.data <- aac.data %>%
  mutate(electionwinner.post = case_when (
    aac.data$electionwinnerpost == 1 ~ 1, 
    aac.data$electionwinnerpost == 2 ~ 0, 
    aac.data$electionwinnerpost == 0 ~ 0))

table(aac.data$electionwinner.post)

##Dominion changed votes (post)## 
aac.data <- aac.data %>%
  mutate(voteschanged.post = case_when (
    aac.data$voteschangedpost == 1 ~ 1, 
    aac.data$voteschangedpost == 2 ~ 0, 
    aac.data$voteschangedpost == 0 ~ 0))

table(aac.data$voteschanged.post)

##Immigrants Voted (post)## 
aac.data <- aac.data %>%
  mutate(illegalvotes.post = case_when (
    aac.data$illegalvotespost == 1 ~ 1, 
    aac.data$illegalvotespost == 2 ~ 0, 
    aac.data$illegalvotespost == 0 ~ 0))

table(aac.data$illegalvotes.post)

##venezuela (post)## 
aac.data <- aac.data %>%
  mutate(venezuela.post= case_when (
    aac.data$venezuelapost == 0 ~ 1, 
    aac.data$venezuelapost == 2 ~ 0, 
    aac.data$venezuelapost == 1 ~ 0))

table(aac.data$venezuela.post)

##uvalde (post)## 
aac.data <- aac.data %>%
  mutate(uvalde.post= case_when (
    aac.data$uvaldepost == 1 ~ 1, 
    aac.data$uvaldepost == 2 ~ 0, 
    aac.data$uvaldepost == 0 ~ 0))

##sandy hook (post)## 
aac.data <- aac.data %>%
  mutate(sandyhook.post= case_when (
    aac.data$sandyhookpost == 0 ~ 0, 
    aac.data$sandyhookpost == 2 ~ 0, 
    aac.data$sandyhookpost == 1 ~ 1))

##Big lie factor variable##

cfamodel2  <- ' h =~ NA*electionwinner.pre + voteschanged.pre + illegalvotes.pre + venezuela.pre
h ~~ 1*h ' 

##Confirmatory factor analysis--model fit indices 
##Results in Appendix B.1##
fit.cfa2 <- sem(cfamodel2, data = aac.data)
summary(fit.cfa2, fit.measures = TRUE, standardized = TRUE)

##Creating New Omnibus Variables##

aac.data$big.lie.pre <- (aac.data$electionwinner.pre + aac.data$voteschanged.pre + aac.data$illegalvotes.pre + 
                           aac.data$venezuela.pre)/4

aac.data$big.lie.post <- (aac.data$electionwinner.post + aac.data$voteschanged.post + aac.data$illegalvotes.post + 
                            aac.data$venezuela.post)/4

##H2: Faith in election security
#2020# 

aac.data$vote2020.scale <- aac.data$vote2020_1 / 100 
table(aac.data$vote2020.scale)

#2024#
aac.data$vote2024.scale <- aac.data$vote2024_1 / 100 
table(aac.data$vote2024.scale)

##H3: Trust in Fox News and Alex Jones

#Fox trust 
aac.data$trust.fox <- aac.data$trustfox_1 / 100 
table(aac.data$trust.fox)

#Jones trust
aac.data$trust.jones <- aac.data$jonestrust_1 / 100 
table(aac.data$trust.jones)

##H4: Future consumption of Fox News and Alex Jones 
##Consumption of Fox News 
aac.data <- aac.data %>%
  mutate(fox.future= case_when (
    aac.data$futuremediachoice_1 == 1 ~ 1, 
    aac.data$futuremediachoice_1 == 2 ~ .5, 
    aac.data$futuremediachoice_1 == 3 ~ 0))

##Consumption of Alex Jones 
aac.data <- aac.data %>%
  mutate(jones.future= case_when (
    aac.data$futurechoicejones_4 == 1 ~ 1, 
    aac.data$futurechoicejones_4 == 2 ~ .5, 
    aac.data$futurechoicejones_4 == 3 ~ 0))

##H5: Future consumption of alternative media sources 

#Consumption of newsmax
aac.data <- aac.data %>%
  mutate(newsmax.future= case_when (
    aac.data$futuremediachoice_2 == 1 ~ 1, 
    aac.data$futuremediachoice_2 == 2 ~ .5, 
    aac.data$futuremediachoice_2 == 3 ~ 0))

#Consumption of OANN 
aac.data <- aac.data %>%
  mutate(oann.future= case_when (
    aac.data$futuremediachoice_3 == 1 ~ 1, 
    aac.data$futuremediachoice_3 == 2 ~ .5, 
    aac.data$futuremediachoice_3 == 3 ~ 0))

#Combined alt-right variable 

aac.data$altright.future <- (aac.data$newsmax.future + aac.data$oann.future)/2

#Consumption of ABC
aac.data <- aac.data %>%
  mutate(abc.future= case_when (
    aac.data$futuremediachoice_5 == 1 ~ 1, 
    aac.data$futuremediachoice_5 == 2 ~ .5, 
    aac.data$futuremediachoice_5 == 3 ~ 0))

#Consumption of CBS 
aac.data <- aac.data %>%
  mutate(cbs.future= case_when (
    aac.data$futuremediachoice_6 == 1 ~ 1, 
    aac.data$futuremediachoice_6 == 2 ~ .5, 
    aac.data$futuremediachoice_6 == 3 ~ 0))

#Combined traditional variable 

aac.data$traditional.future <- (aac.data$abc.future + aac.data$cbs.future)/2

####Hypothesis Testing####
#Note: pre-registration specified that I would conduct simple linear reg
#and regressions w/Lin covariates; models with Lin covariates were main models.
#To restrict replication to only figs/tables that appear in manuscript;
#I show only models w/Lin covariates here. Simple lin models are substantively identical. 

##H1a: admissions reduce support for misperceptions surrounding 2020 election## 

##Regression Model 1: Treatment on Support for Factor/Individual misperceptions 
##w/ Lin adjusted covariates (Appendix B.2)##

models1 <- list(
  "Big Lie"     = lm_lin(big.lie.post ~ treatment.fox, 
                         ~ white + hispanic + college + male + age + fox
                         + media.trust + pid.scale + big.lie.pre,
                         data = aac.data),
  "Trump Won" = lm_lin(electionwinner.post ~ treatment.fox, 
                       ~ white + hispanic + college + male + age + fox
                       + media.trust + pid.scale + electionwinner.pre,
                       data = aac.data),
  "Dominion Changed" = lm_lin(voteschanged.post ~ treatment.fox, 
                              ~ white + hispanic + college + male + age + fox
                              + media.trust + pid.scale + voteschanged.pre,
                              data = aac.data),
  "Immigrants Voted" = lm_lin(illegalvotes.post ~ treatment.fox, 
                              ~ white + hispanic + college + male + age + fox
                              + media.trust + pid.scale + illegalvotes.pre,
                              data = aac.data), 
  "Venezuela Owns" = lm_lin(venezuela.post ~ treatment.fox, 
                            ~ white + hispanic + college + male + age + fox
                            + media.trust + pid.scale + venezuela.pre,
                            data = aac.data))

modelsummary(models1,
             stars = TRUE,
             coef_rename = c("treatment.fox1" = "Admissions",
                             "treatment.fox2" = "Fact check",
                             "white_c" = "white", 
                             "hispanic_c" = "hispanic",
                             "college_c" = "College educated", "male_c" = "male", 
                             "age_c" = "age", "fox_c" = "Fox viewer",
                             "media.trust_c" = "Media trust", 
                             "pid.scale_c" = "PID Scale", 
                             "venezuela.pre_c" = "pre-Ven",
                             "illegalvotes.pre_c" = "pre-illegal",
                             "voteschanged.pre_c" = "pre-votes",
                             "electionwinner.pre_c" = "pre-elec.winner",
                             "big.lie.pre_c" = "pre-big lie"),
             notes = list ('Robust Standard Errors in parentheses',
                           'Source: CloudResearch.'),
             gof_omit = 'DF|Deviance|R2|Log.Lik.|BIC|Std. Errors',
             output = "latex")

##Visualizing H1a (Figure 2a)##
#Renaming
cm <- c('treatment.fox1' = 'Admissions',
        'treatment.fox2' = 'Fact check')

# Shape mapping
shapes <- c('Admissions' = 17,  # Circle
            'Fact check' = 16)   # Triangle

# Background for the plot
b <- list(geom_vline(xintercept = 0, color = 'red'))

# Color specification 
color_list <- c("#000080", "#0000CC", "darkred", "#FF0000", "#004d00", "#009900",
                "purple", "#D8BFD8", "#8B4513", "#CC6600")

b <- list(geom_vline(xintercept = 0, color = 'red'))

# Manually map shapes and colors, ensuring they are applied after renaming
modelplot(models1, coef_map = cm, background = b) +
  aes(color = term, shape = term) +  # Map both color and shape to `term`
  scale_color_manual(values = color_list) +  # Apply custom colors
  scale_shape_manual(values = shapes) +  # Apply custom shapes
  scale_x_continuous(limits = c(-0.20, 0.05)) + 
  theme(
    text = element_text(size = 10),  # Increase text size
    panel.border = element_rect(color = "black", fill = NA, size = .5),  # Add border
    axis.title = element_text(face = "bold"),  # Bold axis titles
    plot.title = element_text(face = "bold")  # Bold plot title
  ) +
  guides(color = guide_legend(reverse = TRUE), 
         shape = guide_legend(reverse = TRUE)) + 
  labs(
    x = 'Change in Endorsement of False Statement', 
    y = '',
    title = 'Effect of Treatment on Election Misperceptions',
    caption = "CloudResearch Connect. Oct. 11-21st, 2024. N = 2382. 
    \n 95% confidence intervals presented with robust standard errors."
  )

##H1b: admissions reduce support for misperceptions surrounding Sandy Hook## 

##Regression Model 2: Treatment on Support for shooting misperceptions 
##w/ Lin adjusted covariates (Appendix B.3)##

models2 <- list(
  "Sandy Hook"     = lm_lin(sandyhook.post ~ treatment.jones, 
                            ~ white + hispanic + college + male + age + jones
                            + media.trust + pid.scale + sandyhook.pre,
                            data = aac.data),
  "Uvalde" = lm_lin(uvalde.post ~ treatment.jones, 
                    ~ white + hispanic + college + male + age + jones
                    + media.trust + pid.scale + uvalde.pre,
                    data = aac.data))

modelsummary(models2,
             stars = TRUE,
             coef_rename = c("treatment.jones1" = "Admissions",
                             "treatment.jones2" = "Fact check",
                             "white_c" = "white",
                             "hispanic_c" = "hispanic",
                             "college_c" = "College educated", "male_c" = "male", 
                             "age_c" = "age", "jones_c" = "Jones viewer",
                             "media.trust_c" = "Media trust", "pid.scale_c" = "PID Scale",
                             "sandyhook.pre_c" = "pre-Sandy Hook",
                             "uvalde.pre_c" = "pre_Uvalde"),
             notes = list ('Robust Standard Errors in parentheses',
                           'Source: CloudResearch.'),
             gof_omit = 'DF|Deviance|R2|Log.Lik.|BIC|Std. Errors',
             output = "latex")

##Visualizing H1b (Fig 2b)##
cm <- c('treatment.jones1' = 'Admissions',
        'treatment.jones2' = 'Fact check')

# Shape mapping
shapes <- c('Admissions' = 17,  # Circle
            'Fact check' = 16)   # Triangle

# Background for the plot
b <- list(geom_vline(xintercept = 0, color = 'red'))

# Color specification 
color_list <- c("#FF0000", "#004d00")

b <- list(geom_vline(xintercept = 0, color = 'red'))

# Manually map shapes and colors, ensuring they are applied after renaming
modelplot(models2, coef_map = cm, background = b) +
  aes(color = term, shape = term) +  # Map both color and shape to `term`
  scale_color_manual(values = color_list) +  # Apply custom colors
  scale_shape_manual(values = shapes) +  # Apply custom shapes
  scale_x_continuous(limits = c(-0.20, 0.05)) + #Apply custom X axis 
  theme(
    text = element_text(size = 10),  # Increase text size
    panel.border = element_rect(color = "black", fill = NA, size = .5),  # Add border
    axis.title = element_text(face = "bold"),  # Bold axis titles
    plot.title = element_text(face = "bold")  # Bold plot title
  ) +
  guides(color = guide_legend(reverse = TRUE), shape = guide_legend(reverse = TRUE)) + 
  labs(
    x = 'Change in Endorsement of False Statement', 
    y = '',
    title = 'Effect of Treatment on Shooting Misperceptions',
    caption = "CloudResearch Connect. Oct. 11-21st, 2024. N = 2382. 
    \n 95% confidence intervals presented with robust standard errors."
  )

#H2a: treatments increase faith in election security 

##Regression Model 3: Treatment on election faith 
##w/ Lin adjusted covariates (Appendix B.4)##

models3 <- list(
  "Election 2020"     = lm_lin(vote2020.scale ~ treatment.fox, 
                               ~ white + hispanic + college + male + age + fox
                               + media.trust + pid.scale,
                               data = aac.data),
  "Election 2024" = lm_lin(vote2024.scale ~ treatment.fox, 
                           ~ white + hispanic + college + male + age + fox
                           + media.trust + pid.scale,
                           data = aac.data))

modelsummary(models3,
             stars = TRUE,
             coef_rename = c("treatment.fox1" = "Admissions",
                             "treatment.fox2" = "Fact check",
                             "white_c" = "white", 
                             "hispanic_c" = "hispanic",
                             "college_c" = "College educated", "male_c" = "male", 
                             "age_c" = "age", "fox_c" = "Fox viewer",
                             "media.trust_c" = "Media trust", "pid.scale_c" = "PID Scale"),
             notes = list ('Robust Standard Errors in parentheses',
                           'Source: CloudResearch.'),
             gof_omit = 'DF|Deviance|R2|Log.Lik.|BIC|Std. Errors',
             output = "latex")

#Visualizing H3a (Figure 3a)

cm <- c('treatment.fox1' = 'Admissions',
        'treatment.fox2' = 'Fact check')

# Shape mapping
shapes <- c('Admissions' = 17,  # Circle
            'Fact check' = 16)   # Triangle

# Background for the plot
b <- list(geom_vline(xintercept = 0, color = 'red'))

# Color specification 
color_list <- c("#FF0000", "#004d00")

b <- list(geom_vline(xintercept = 0, color = 'red'))

# Manually map shapes and colors, ensuring they are applied after renaming
modelplot(models3, coef_map = cm, background = b) +
  aes(color = term, shape = term) +  # Map both color and shape to `term`
  scale_color_manual(values = color_list) +  # Apply custom colors
  scale_shape_manual(values = shapes) +  # Apply custom shapes
  theme(
    text = element_text(size = 10),  # Increase text size
    panel.border = element_rect(color = "black", fill = NA, size = .5),  # Add border
    axis.title = element_text(face = "bold"),  # Bold axis titles
    plot.title = element_text(face = "bold")  # Bold plot title
  ) +
  guides(color = guide_legend(reverse = TRUE), shape = guide_legend(reverse = TRUE)) + 
  labs(
    x = 'Change in Faith in Election Security', 
    y = '',
    title = 'Effect of Treatment on Election Trust',
    caption = "CloudResearch Connect. Oct. 11-21st, 2024. N = 2336-2372. \n 95% confidence intervals presented with robust standard errors."
  )

#H3a and H3b: admissions will decrease faith in Fox Nes and Alex Jones 

##Regression Model 4: effect of admissions on trust in Fox News and Alex Jones  
##w/ Lin adjusted covariates (Appendix B.5)##

models4 <- list(
  "Fox Trust"     = lm_lin(trust.fox ~ treatment.fox, 
                           ~ white + hispanic + college + male + age + fox
                           + media.trust + pid.scale,
                           data = aac.data),
  "Jones Trust" = lm_lin(trust.jones ~ treatment.jones, 
                         ~ white + hispanic + college + male + age + jones
                         + media.trust + pid.scale,
                         data = aac.data))

modelsummary(models4,
             stars = TRUE,
             coef_rename = c("treatment.fox1" = "Admissions (Fox)",
                             "treatment.fox2" = "Fact check (Fox)",
                             "treatment.jones1" = "Admissions (Jones)",
                             "treatment.jones2" = "Fact check (Jones)",
                             "white_c" = "white", 
                             "hispanic_c" = "hispanic",
                             "college_c" = "College educated", "male_c" = "male", 
                             "age_c" = "age", "fox_c" = "Fox viewer",
                             "jones_c" = "Jones viewer",
                             "media.trust_c" = "Media trust", "pid.scale_c" = "PID Scale"),
             notes = list ('Robust Standard Errors in parentheses',
                           'Source: CloudResearch.'),
             gof_omit = 'DF|Deviance|R2|Log.Lik.|BIC|Std. Errors',
             output = "latex")

#H4a and H4b: admissions will decrease intended consumption of Fox News 
## and Alex Jones w/ Lin adjusted covariates (Appendix B.6)

##Regression model 5: effect of admissions on future consumption of
#Fox and Jones (Lin covariates)

models5 <- list(
  "Future Fox"     = lm_lin(fox.future ~ treatment.fox, 
                            ~ white + hispanic + college + male + age + fox
                            + media.trust + pid.scale,
                            data = aac.data),
  "Future Jones" = lm_lin(jones.future ~ treatment.jones, 
                          ~ white + hispanic + college + male + age + jones
                          + media.trust + pid.scale,
                          data = aac.data))

modelsummary(models5,
             stars = TRUE,
             coef_rename = c("treatment.fox1" = "Admissions (Fox)",
                             "treatment.fox2" = "Fact check (Fox)",
                             "treatment.jones1" = "Admissions (Jones)",
                             "treatment.jones2" = "Fact check (Jones)",
                             "white_c" = "white", 
                             "hispanic_c" = "hispanic",
                             "college_c" = "College educated", "male_c" = "male", 
                             "age_c" = "age", "fox_c" = "Fox viewer",
                             "jones_c" = "Jones viewer",
                             "media.trust_c" = "Media trust", "pid.scale_c" = "PID Scale"),
             notes = list ('Robust Standard Errors in parentheses',
                           'Source: CloudResearch.'),
             gof_omit = 'DF|Deviance|R2|Log.Lik.|BIC|Std. Errors',
             output = "latex")

##Visualizing H3 & H4 for output (Figure 3b)

models6 <- list(
  "Future Fox"     = lm_lin(fox.future ~ treatment.fox, 
                            ~ white + hispanic + college + male + age + fox
                            + media.trust + pid.scale,
                            data = aac.data),
  "Future Jones" = lm_lin(jones.future ~ treatment.jones, 
                          ~ white + hispanic + college + male + age + jones
                          + media.trust + pid.scale,
                          data = aac.data),
  "Fox Trust"     = lm_lin(trust.fox ~ treatment.fox, 
                           ~ white + hispanic + college + male + age + fox
                           + media.trust + pid.scale,
                           data = aac.data),
  "Jones Trust" = lm_lin(trust.jones ~ treatment.jones, 
                         ~ white + hispanic + college + male + age + jones
                         + media.trust + pid.scale,
                         data = aac.data))

cm <- c('treatment.fox1' = 'Admissions',
        'treatment.fox2' = 'Fact check',
        'treatment.jones1' = 'Admissions',
        'treatment.jones2' = 'Fact check')

# Shape mapping
shapes <- c('Admissions' = 17,  # Circle
            'Fact check' = 16)   # Triangle

# Background for the plot
b <- list(geom_vline(xintercept = 0, color = 'red'))

# Color specification 
color_list <- c("#FF0000", "#004d00", "#0000CC","purple") 

b <- list(geom_vline(xintercept = 0, color = 'red'))

# Manually map shapes and colors, ensuring they are applied after renaming
modelplot(models6, coef_map = cm, background = b) +
  aes(color = term, shape = term) +  # Map both color and shape to `term`
  scale_color_manual(values = color_list) +  # Apply custom colors
  scale_shape_manual(values = shapes) +  # Apply custom shapes
  theme(
    text = element_text(size = 10),  # Increase text size
    panel.border = element_rect(color = "black", fill = NA, size = .5),  # Add border
    axis.title = element_text(face = "bold"),  # Bold axis titles
    plot.title = element_text(face = "bold")  # Bold plot title
  ) +
  guides(color = guide_legend(reverse = TRUE), shape = guide_legend(reverse = TRUE)) + 
  labs(
    x = 'Change in Trust and Future Consumption of Fox/Jones', 
    y = '',
    title = 'Effect of Treatment on Media DVs',
    caption = "CloudResearch Connect. Oct. 11-21st, 2024. N = 2370-2379. \n 95% confidence intervals presented with robust standard errors."
  )


#H5a and h6a: treatments affect consumption of alt-right/traditional media 

##Regression Model 7: Treatment on intended consumption of alt-right/trad media 
##w/ Lin adjusted covariates (Appendix B.7)##

models7 <- list(
  "Future Alt Right"     = lm_lin(altright.future ~ treatment.fox, 
                                  ~ white + hispanic + college + male + age + fox
                                  + media.trust + pid.scale,
                                  data = aac.data),
  "Future Traditional" = lm_lin(traditional.future ~ treatment.fox, 
                                ~ white + hispanic + college + male + age + fox
                                + media.trust + pid.scale,
                                data = aac.data))

modelsummary(models7,
             stars = TRUE,
             coef_rename = c("treatment.fox1" = "Admissions",
                             "treatment.fox2" = "Fact check",
                             "white_c" = "white", 
                             "hispanic_c" = "hispanic",
                             "college_c" = "College educated", "male_c" = "male", 
                             "age_c" = "age", "fox_c" = "Fox viewer",
                             "media.trust_c" = "Media trust", "pid.scale_c" = "PID Scale"),
             notes = list ('Robust Standard Errors in parentheses',
                           'Source: CloudResearch.'),
             gof_omit = 'DF|Deviance|R2|Log.Lik.|BIC|Std. Errors',
             output = "latex")

##Two one-sided equivalence test--are the effects negligible?##
##Note: included in original (Stage 1) submission as part of JOP mandatory 
##design table for all non-significant results, included here for completeness

##Testing whether means are precisely null
##Calculate mean and SD of altright consumption 
aac.data %>%
  group_by(treatment.fox) %>%
  summarise(
    mean_altright = mean(altright.future, na.rm = TRUE),
    sd_altright = sd(altright.future, na.rm = TRUE),
    n = n()
  )

##Using TOSTer to calculate equivalence (admissions to control)##
TOSTtwo(m1=.175, m2=.184, sd1=0.284, 
        sd2=0.288, n1=828, n2=796, 
        low_eqbound_d=-0.05, high_eqbound_d=0.05, 
        alpha = 0.05, var.equal = TRUE)

##Using TOSTer to calculate equivalence (fact-check to control)
TOSTtwo(m1=.175, m2=.192, sd1=0.28, 
        sd2=0.286, n1=828, n2=760, 
        low_eqbound_d=-0.05, high_eqbound_d=0.05, 
        alpha = 0.05, var.equal = TRUE)

##Findings--effect is statistically equivalent to zero for both treatment

##Calculate mean and SD of traditional future consumption##
aac.data %>%
  group_by(treatment.fox) %>%
  summarise(
    mean_altright = mean(traditional.future, na.rm = TRUE),
    sd_altright = sd(traditional.future, na.rm = TRUE),
    n = n()
  )

##Using TOSTer to calculate equivalence (admissions to control)##
TOSTtwo(m1=.301, m2=.299, sd1=0.351, 
        sd2=0.348, n1=828, n2=796, 
        low_eqbound_d=-0.05, high_eqbound_d=0.05, 
        alpha = 0.05, var.equal = TRUE)

##Using TOSTer to calculate equivalence (fact-check to control)
TOSTtwo(m1=.301, m2=.299, sd1=0.351, 
        sd2=0.347, n1=828, n2=760, 
        low_eqbound_d=-0.05, high_eqbound_d=0.05, 
        alpha = 0.05, var.equal = TRUE)

##Finding--effect is statistically equivalent to zero for both treatments

#H7a: effectiveness of admissions on election misperceptions w/fact-check as
#baseline

##First set fact check as baseline 

##Experiment on Fox News## 
aac.data <- aac.data %>%
  mutate(treatment.fox.fact.baseline = case_when (
    aac.data$FL_22_DO == "FL_23" ~ 2, #control
    aac.data$FL_22_DO == "FL_24" ~ 1, #admissions
    aac.data$FL_22_DO == "FL_63" ~ 0)) #fact check

###Regression Model 8: Treatment on Support for election misperceptions w/fact
##check as baseline and Lin adjusted covariates (Appendix B.8)##

models8 <- list(
  "Big Lie"     = lm_lin(big.lie.post ~ treatment.fox.fact.baseline, 
                         ~ white + hispanic + college + male + age + fox
                         + media.trust + pid.scale + big.lie.pre,
                         data = aac.data),
  "Trump Won" = lm_lin(electionwinner.post ~ treatment.fox.fact.baseline, 
                       ~ white + hispanic + college + male + age + fox
                       + media.trust + pid.scale + electionwinner.pre,
                       data = aac.data),
  "Dominion Changed" = lm_lin(voteschanged.post ~ treatment.fox.fact.baseline, 
                              ~ white + hispanic + college + male + age + fox
                              + media.trust + pid.scale + voteschanged.pre,
                              data = aac.data),
  "Immigrants Voted" = lm_lin(illegalvotes.post ~ treatment.fox.fact.baseline, 
                              ~ white + hispanic + college + male + age + fox
                              + media.trust + pid.scale + illegalvotes.pre,
                              data = aac.data), 
  "Venezuela Owns" = lm_lin(venezuela.post ~ treatment.fox.fact.baseline, 
                            ~ white + hispanic +  college + male + age + fox
                            + media.trust + pid.scale + venezuela.pre,
                            data = aac.data))

modelsummary(models8,
             stars = TRUE,
             coef_rename = c("treatment.fox.fact.baseline1" = "Admissions",
                             "treatment.fox.fact.baseline2" = "Control",
                             "white_c" = "white", 
                             "hispanic_c" = "hispanic",
                             "college_c" = "College educated", "male_c" = "male", 
                             "age_c" = "age", "fox_c" = "Fox viewer",
                             "media.trust_c" = "Media trust", "pid.scale_c" = "PID Scale",
                             "venezuela.pre_c" = "pre-Ven",
                             "illegalvotes.pre_c" = "pre-illegal",
                             "voteschanged.pre_c" = "pre-votes",
                             "electionwinner.pre_c" = "pre-elec.winner",
                             "big.lie.pre_c" = "pre-big lie"),
             notes = list ('Robust Standard Errors in parentheses',
                           'Source: CloudResearch.'),
             gof_omit = 'DF|Deviance|R2|Log.Lik.|BIC|Std. Errors',
             output = "latex")

##H7b: testing effectiveness of treatment on misperceptions w/fact check baseline 
#(Sandy Hook)

#First set fact-check as baseline
##Experiment on Sandy Hook## 
aac.data <- aac.data %>%
  mutate(treatment.jones.fact.baseline = case_when (
    aac.data$FL_61_DO == "FL_62" ~ 2, #control
    aac.data$FL_61_DO == "FL_64" ~ 1, #admissions
    aac.data$FL_61_DO == "FL_65" ~ 0)) #fact check

#Regression Model 9: Treatment on Support for shooting misperceptions 
#w/Lin covariates (Appendix B.9)

models9 <- list(
  "Sandy Hook"     = lm_lin(sandyhook.post ~ treatment.jones.fact.baseline, 
                            ~ white + hispanic + college + male + age + jones
                            + media.trust + pid.scale + sandyhook.pre,
                            data = aac.data),
  "Uvalde" = lm_lin(uvalde.post ~ treatment.jones.fact.baseline, 
                    ~ white + college + hispanic + male + age + jones
                    + media.trust + pid.scale + uvalde.pre,
                    data = aac.data))

modelsummary(models9,
             stars = TRUE,
             coef_rename = c("treatment.jones.fact.baseline1" = "Admissions",
                             "treatment.jones.fact.baseline2" = "Control",
                             "white_c" = "white", 
                             "hispanic_c" = "hispanic",
                             "college_c" = "College educated", "male_c" = "male", 
                             "age_c" = "age", "jones_c" = "Jones viewer",
                             "media.trust_c" = "Media trust", "pid.scale_c" = "PID Scale",
                             "sandyhook.pre_c" = "pre-Sandy Hook",
                             "uvalde.pre_c" = "pre-Uvalde"),
             notes = list ('Robust Standard Errors in parentheses',
                           'Source: CloudResearch.'),
             gof_omit = 'DF|Deviance|R2|Log.Lik.|BIC|Std. Errors',
             output = "latex")

####Differential attrition test####
##NOTE: code will not work if you run lines 95-119, which drops all non-completes##
##To run, skip those lines## 

##Completion rates by experimental treatment## 

compute_completion_rate <- function(data, treatment_var) {
  data %>%
    group_by(!!sym(treatment_var)) %>%
    summarise(Completion_Rate = mean(complete, na.rm = TRUE)) %>%
    mutate(Completion_Rate = round(Completion_Rate * 100, 1)) # Convert to percentage
}

# Compute completion rates
jones_completion <- compute_completion_rate(aac.data, "treatment.jones")
fox_completion <- compute_completion_rate(aac.data, "treatment.fox")

cat("\nCompletion Rate by Treatment.Jones:\n")
print(jones_completion)

cat("\nCompletion Rate by Treatment.Fox:\n")
print(fox_completion)

##Regress completion on treatment condition for both experiments (Appendix B.10)
attrition.models <- list(
  "Attrition Sandy Hook"     = lm_lin(complete ~ treatment.jones, 
                                      ~ white + hispanic + college + male + age 
                                      + jones + media.trust + pid.scale,
                                      data = aac.data),
  "Attrition Dominion" = lm_lin(complete ~ treatment.fox, 
                                ~ white + college + hispanic + male + age + 
                                  fox + media.trust + pid.scale,
                                data = aac.data))

modelsummary(attrition.models,
             stars = TRUE,
             coef_rename = c("treatment.fox1" = "Admissions",
                             "treatment.fox2" = "Fact check",
                             "treatment.jones1" = "Admissions",
                             "treatment.jones2" = "Fact check",
                             "white_c" = "white", 
                             "hispanic_c" = "hispanic",
                             "college_c" = "College educated", 
                             "male_c" = "male", 
                             "age_c" = "age", 
                             "fox_c" = "Fox viewer",
                             "media.trust_c" = "Media trust", 
                             "pid.scale_c" = "PID Scale"),
             notes = list ('Robust Standard Errors in parentheses',
                           'Source: CloudResearch.'),
             gof_omit = 'DF|Deviance|R2|Log.Lik.|BIC|Std. Errors',
             output = "latex")

####Non-prereigstered analyses####

#Calculating IV levels (misperceptions and media) for each group 
#Calculate beliefs by treatment group for misperceptions##
#Used to contextualize effect size (calculating percentage change) 
#in discussion of final manuscript 
compute_misperception_means <- 
  function(data, treatment_var, misperception_vars) {
  data %>%
    group_by(!!sym(treatment_var)) %>%
    summarise(across(all_of(misperception_vars), mean, na.rm = TRUE)) %>%
    mutate(across(where(is.numeric), round, 3))  # Round to 2 decimal places
}

# Misperception variables for each treatment
fox_misperceptions <- c("electionwinner.post", "voteschanged.post", "venezuela.post",  "illegalvotes.post")
jones_misperceptions <- c("uvalde.post", "sandyhook.post")

# Compute means for each treatment
fox_means <- compute_misperception_means(aac.data, "treatment.fox", fox_misperceptions)
jones_means <- compute_misperception_means(aac.data, "treatment.jones", jones_misperceptions)

# Print results
cat("\nMisperception Beliefs by Treatment.Fox:\n")
print(fox_means)

cat("\nMisperception Beliefs by Treatment.Jones:\n")
print(jones_means)

##Trust and consumption media variables##
compute_means <- function(data, treatment_var, variables) {
  data %>%
    group_by(!!sym(treatment_var)) %>%
    summarise(across(all_of(variables), mean, na.rm = TRUE)) %>%
    mutate(across(where(is.numeric), round, 3))  # Round to 2 decimal places
}

# Variables to analyze
fox_vars <- c("trust.fox", "fox.future")
jones_vars <- c("trust.jones", "jones.future")

# Compute means for each treatment
fox_means <- compute_means(aac.data, "treatment.fox", fox_vars)
jones_means <- compute_means(aac.data, "treatment.jones", jones_vars)

# Print results
cat("\nTrust and Future Beliefs by Treatment.Fox:\n")
print(fox_means)

cat("\nTrust and Future Beliefs by Treatment.Jones:\n")
print(jones_means)

#Exploratory analysis: is admissions sig more effective than fact-check 
#at decreasing trust/future consumption of Jones/Fox? (Appendix B.11)

exploratory_hyp <- list(
  "Trust Fox"     = lm_lin(trust.fox ~ treatment.fox.fact.baseline, 
                         ~ white + hispanic + college + male + age + fox
                         + media.trust + pid.scale,
                         data = aac.data),
  "Fox Future" = lm_lin(fox.future ~ treatment.fox.fact.baseline, 
                       ~ white + hispanic + college + male + age + fox
                       + media.trust + pid.scale,
                       data = aac.data),
  "Trust Jones" = lm_lin(trust.jones ~ treatment.jones.fact.baseline, 
                              ~ white + hispanic + college + male + age + jones
                              + media.trust + pid.scale,
                              data = aac.data),
  "Jones Future" = lm_lin(jones.future ~ treatment.jones.fact.baseline, 
                              ~ white + hispanic + college + male + age + jones
                              + media.trust + pid.scale,
                              data = aac.data))

modelsummary(exploratory_hyp,
             stars = TRUE,
             coef_rename = c("treatment.fox.fact.baseline1" = "Admissions",
                             "treatment.fox.fact.baseline2" = "Control",
                             "treatment.jones.fact.baseline1" = "Admissions",
                             "treatment.jones.fact.baseline2" = "Control",
                             "white_c" = "white", 
                             "hispanic_c" = "hispanic",
                             "college_c" = "College educated", "male_c" = "male", 
                             "age_c" = "age", "fox_c" = "Fox viewer",
                             "jones_c" = "Jones viewer",
                             "media.trust_c" = "Media trust", "pid.scale_c" = "PID Scale",
                             "venezuela.pre_c" = "pre-Ven"),
             notes = list ('Robust Standard Errors in parentheses',
                           'Source: CloudResearch.'),
             gof_omit = 'DF|Deviance|R2|Log.Lik.|BIC|Std. Errors',
             output = "latex")



